Macrophages derived from Control (CT) and PD individuals were exposed to myelin prepared from CT individuals.
df = read_xlsx("ctrl_myelin-macrophages.xlsx")
df
df %<>%
pivot_longer(-TIME, values_to = "INT", names_to = "ID") %>%
separate(ID, sep = 2, into = "STATUS", remove = FALSE) %>%
select(ID, STATUS, TIME, INT) %>%
arrange(ID, TIME)
df
ggplot(data = df, mapping = aes(x = TIME, y = INT, color = STATUS)) + geom_smooth(method = "loess", formula="y ~ x")
ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = STATUS)) +
geom_point() +
geom_line()
ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = STATUS)) + geom_smooth(method = "lm", formula="y ~ 0 + x", se = FALSE)
m1.fit = lmer(INT ~ 0 + TIME*STATUS + (0 + TIME | ID), data=df)
Warning in as_lmerModLT(model, devfun) :
Model may not have converged with 1 eigenvalue close to zero: 1.3e-11
summary(m1.fit)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ 0 + TIME * STATUS + (0 + TIME | ID)
Data: df
REML criterion at convergence: 14800
Scaled residuals:
Min 1Q Median 3Q Max
-3.152 -0.443 -0.061 0.439 7.594
Random effects:
Groups Name Variance Std.Dev.
ID TIME 13955019515874 3735642
Residual 120010088567134 10954912
Number of obs: 420, groups: ID, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
TIME 6081194.1 1132610.7 18.3 5.37 0.000039822857 ***
STATUSCT 17763784.6 1391548.9 114154268801901680132096.0 12.77 < 0.0000000000000002 ***
STATUSPD 10266197.3 1538415.2 79120202192516307484672.0 6.67 0.000000000025 ***
TIME:STATUSPD -2356761.2 1688396.3 18.3 -1.40 0.18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
TIME STATUSC STATUSP
STATUSCT -0.090
STATUSPD 0.000 0.000
TIME:STATUS -0.671 0.060 -0.067
m1.res = simulateResiduals(fittedModel = m1.fit)
plot(m1.res)
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
m2.fit = lmer(INT ~ 0 + TIME*STATUS + (0 + TIME | ID), data=df %>% filter(TIME <= 10))
Warning in as_lmerModLT(model, devfun) :
Model may not have converged with 1 eigenvalue close to zero: 1.2e-11
summary(m2.fit)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ 0 + TIME * STATUS + (0 + TIME | ID)
Data: df %>% filter(TIME <= 10)
REML criterion at convergence: 7589
Scaled residuals:
Min 1Q Median 3Q Max
-2.268 -0.413 -0.039 0.234 10.734
Random effects:
Groups Name Variance Std.Dev.
ID TIME 18017163807632 4244663
Residual 64442444229858 8027605
Number of obs: 220, groups: ID, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
TIME 8716594.4 1300454.6 18.8 6.70 0.0000022 ***
STATUSCT 5488175.9 1365297.6 1883035153459291881472.0 4.02 0.0000583 ***
STATUSPD 4298811.9 1509393.3 420829932100792145149952.0 2.85 0.0044 **
TIME:STATUSPD -3747708.4 1938603.3 18.8 -1.93 0.0684 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
TIME STATUSC STATUSP
STATUSCT -0.150
STATUSPD 0.000 0.000
TIME:STATUS -0.671 0.101 -0.111
m2.res = simulateResiduals(fittedModel = m2.fit)
plot(m2.res)
ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = STATUS)) + geom_smooth(method = "lm", formula="y ~ x", se = FALSE)
m3.fit = lmer(INT ~ TIME*STATUS + (TIME | ID), data=df)
Warning in as_lmerModLT(model, devfun) :
Model may not have converged with 1 eigenvalue close to zero: 1.4e-11
summary(m3.fit)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ TIME * STATUS + (TIME | ID)
Data: df
REML criterion at convergence: 14777
Scaled residuals:
Min 1Q Median 3Q Max
-3.257 -0.338 -0.020 0.384 7.875
Random effects:
Groups Name Variance Std.Dev. Corr
ID (Intercept) 46033779618314 6784820
TIME 12199710583965 3492808 0.44
Residual 108281639341729 10405846
Number of obs: 420, groups: ID, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17763785 2435581 18 7.29 0.00000089 ***
TIME 6081194 1059174 18 5.74 0.00001922 ***
STATUSPD -7497587 3630750 18 -2.07 0.054 .
TIME:STATUSPD -2356761 1578923 18 -1.49 0.153
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) TIME STATUS
TIME 0.315
STATUSPD -0.671 -0.211
TIME:STATUS -0.211 -0.671 0.315
m3.res = simulateResiduals(fittedModel = m3.fit)
plot(m3.res)
m4.fit = lmer(INT ~ TIME*STATUS + (TIME | ID), data=df %>% filter(TIME <= 10))
Warning in as_lmerModLT(model, devfun) :
Model may not have converged with 1 eigenvalue close to zero: 1.2e-11
summary(m4.fit)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ TIME * STATUS + (TIME | ID)
Data: df %>% filter(TIME <= 10)
REML criterion at convergence: 7587
Scaled residuals:
Min 1Q Median 3Q Max
-2.642 -0.387 -0.020 0.213 10.810
Random effects:
Groups Name Variance Std.Dev. Corr
ID (Intercept) 9830589540884 3135377
TIME 17362487055394 4166832 0.12
Residual 61633698670715 7850713
Number of obs: 220, groups: ID, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5488176 1635996 18 3.35 0.0035 **
TIME 8716594 1276458 18 6.83 0.0000022 ***
STATUSPD -1189364 2438799 18 -0.49 0.6317
TIME:STATUSPD -3747708 1902831 18 -1.97 0.0645 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) TIME STATUS
TIME -0.052
STATUSPD -0.671 0.035
TIME:STATUS 0.035 -0.671 -0.052
m4.res = simulateResiduals(fittedModel = m4.fit)
plot(m4.res)
ggplot(data = df, mapping = aes(x = factor(TIME), y = INT, color = STATUS)) +
geom_boxplot()
summary(aov(INT ~ STATUS*factor(TIME) + Error(ID), data = df))
Error: ID
Df Sum Sq Mean Sq F value Pr(>F)
STATUS 1 100316595603111792 100316595603111792 3.23 0.089 .
Residuals 18 558649433744305344 31036079652461408
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
factor(TIME) 20 405726546818305024 20286327340915252 37.89 <0.0000000000000002 ***
STATUS:factor(TIME) 20 23064109848552772 1153205492427638 2.15 0.003 **
Residuals 360 192757182833994944 535436618983319
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(lmer(INT ~ STATUS*factor(TIME) + (1 | ID), data=df))
Warning in as_lmerModLT(model, devfun) :
Model may not have converged with 1 eigenvalue close to zero: 2.7e-12
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
STATUS 1729285478721196 1729285478721196 1 18 3.23 0.089 .
factor(TIME) 382825994343270720 19141299717163536 20 Inf 35.75 <0.0000000000000002 ***
STATUS:factor(TIME) 23064109848552360 1153205492427618 20 Inf 2.15 0.002 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(lmer(INT ~ STATUS*factor(TIME) + (1 | ID), data=df))
Warning in as_lmerModLT(model, devfun) :
Model may not have converged with 1 eigenvalue close to zero: 2.7e-12
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: INT ~ STATUS * factor(TIME) + (1 | ID)
Data: df
REML criterion at convergence: 14062
Scaled residuals:
Min 1Q Median 3Q Max
-2.528 -0.571 -0.004 0.487 3.333
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 1452411626813050 38110519
Residual 535436618009765 23139503
Number of obs: 420, groups: ID, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 31620.1 13442971.3 32.6 0.00 0.99814
STATUSPD -21728.0 20039598.4 32.6 0.00 0.99914
factor(TIME)1 8781401.4 9866717.4 42304711666901590016.0 0.89 0.37346
factor(TIME)2 23768179.7 9866717.4 41019441770847150080.0 2.41 0.01600 *
factor(TIME)3 34487435.5 9866717.4 43569307823568625664.0 3.50 0.00047 ***
factor(TIME)4 50513832.0 9866717.4 38613191227455332352.0 5.12 0.00000030615385 ***
factor(TIME)5 51811985.3 9866717.4 45441836676151836672.0 5.25 0.00000015112156 ***
factor(TIME)6 58989519.0 9866717.4 40347422911320178688.0 5.98 0.00000000225013 ***
factor(TIME)7 66618580.2 9866717.4 40132082484553932800.0 6.75 0.00000000001460 ***
factor(TIME)8 73883423.8 9866717.4 38849758612996620288.0 7.49 0.00000000000007 ***
factor(TIME)9 82034966.7 9866717.4 39666229042097963008.0 8.31 < 0.0000000000000002 ***
factor(TIME)10 88545483.7 9866717.4 38580371851014479872.0 8.97 < 0.0000000000000002 ***
factor(TIME)11 92236016.1 9866717.4 43039447898024665088.0 9.35 < 0.0000000000000002 ***
factor(TIME)12 95729228.8 9866717.4 38863621046500278272.0 9.70 < 0.0000000000000002 ***
factor(TIME)13 102929904.9 9866717.4 40508903641680633856.0 10.43 < 0.0000000000000002 ***
factor(TIME)14 106728537.7 9866717.4 41475839895920320512.0 10.82 < 0.0000000000000002 ***
factor(TIME)15 111739025.5 9866717.4 43975724117093941248.0 11.32 < 0.0000000000000002 ***
factor(TIME)16 112668644.3 9866717.4 41299628101740830720.0 11.42 < 0.0000000000000002 ***
factor(TIME)17 117035519.1 9866717.4 42709417196245622784.0 11.86 < 0.0000000000000002 ***
factor(TIME)18 119225475.0 9866717.4 41683750762701217792.0 12.08 < 0.0000000000000002 ***
factor(TIME)19 123280234.6 9866717.4 40807518497940914176.0 12.49 < 0.0000000000000002 ***
factor(TIME)20 128418824.0 9866717.4 42496879823159566336.0 13.02 < 0.0000000000000002 ***
STATUSPD:factor(TIME)1 -2192325.3 14708433.9 302803913031962787840.0 -0.15 0.88151
STATUSPD:factor(TIME)2 -7940390.1 14708433.9 263360823376362176512.0 -0.54 0.58930
STATUSPD:factor(TIME)3 -11572207.9 14708433.9 315759423867243003904.0 -0.79 0.43141
STATUSPD:factor(TIME)4 -23749022.1 14708433.9 286506603967034687488.0 -1.61 0.10639
STATUSPD:factor(TIME)5 -20600772.6 14708433.9 311933493985617313792.0 -1.40 0.16133
STATUSPD:factor(TIME)6 -23914955.2 14708433.9 321361150475973558272.0 -1.63 0.10396
STATUSPD:factor(TIME)7 -27324733.9 14708433.9 243030239043486711808.0 -1.86 0.06320 .
STATUSPD:factor(TIME)8 -30697005.5 14708433.9 242991815025575854080.0 -2.09 0.03689 *
STATUSPD:factor(TIME)9 -33806339.4 14708433.9 279563481426973523968.0 -2.30 0.02154 *
STATUSPD:factor(TIME)10 -37170207.8 14708433.9 241790702101374074880.0 -2.53 0.01150 *
STATUSPD:factor(TIME)11 -39111121.5 14708433.9 289315740243068125184.0 -2.66 0.00784 **
STATUSPD:factor(TIME)12 -37386024.8 14708433.9 235730788509275062272.0 -2.54 0.01103 *
STATUSPD:factor(TIME)13 -40031013.0 14708433.9 242992509354682646528.0 -2.72 0.00650 **
STATUSPD:factor(TIME)14 -38876187.6 14708433.9 270572760748935938048.0 -2.64 0.00821 **
STATUSPD:factor(TIME)15 -42234173.8 14708433.9 313696440606546329600.0 -2.87 0.00409 **
STATUSPD:factor(TIME)16 -43377542.1 14708433.9 273472501187227844608.0 -2.95 0.00319 **
STATUSPD:factor(TIME)17 -46923491.4 14708433.9 306408864699436498944.0 -3.19 0.00142 **
STATUSPD:factor(TIME)18 -47063593.3 14708433.9 292368502840325636096.0 -3.20 0.00138 **
STATUSPD:factor(TIME)19 -48010183.6 14708433.9 271948042683244412928.0 -3.26 0.00110 **
STATUSPD:factor(TIME)20 -49931609.0 14708433.9 284229729013514764288.0 -3.39 0.00069 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation matrix not shown by default, as p = 42 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
m = df %>%
nest(data = -c(ID, STATUS)) %>%
mutate(
model = map(.x = data, ~ drm(.x$INT ~ .x$TIME, fct = DRC.negExp())),
tidied = map(model, tidy)
)
m
tt = m %>%
unnest(tidied) %>%
pivot_wider(id_cols = c(ID, STATUS), names_from = term, values_from = estimate) %>%
rename(PLATEAU = a, RATE = c)
tt
t.test(PLATEAU ~ STATUS, data = tt)
Welch Two Sample t-test
data: PLATEAU by STATUS
t = 0.39, df = 15, p-value = 0.7
alternative hypothesis: true difference in means between group CT and group PD is not equal to 0
95 percent confidence interval:
-130349196 188986685
sample estimates:
mean in group CT mean in group PD
175235687 145916942
ggplot(data = df, mapping = aes(x = TIME, y = INT, group = ID, color = ID)) +
geom_point() +
geom_line()
See also: